New special functions solving nonlinear 
autonomous dynamical systems 



Leon Brenig0 

Abstract 

A general solution is found for a large class of time continuous au- 
tonomous nonlinear dynamical systems, the so-called quasi-polynomial 
systems. This solution is expressed in terms of a new type of spe- 
cial functions defined via their Taylor series. The coefficients of these 
Taylor series are generated by a tensor that generalizes the factorial 
function and has a combinatorial meaning. The existence of these func- 
tions raises the question of the relation between them and the chaotic 
behaviour of the solutions that may appear for the quasi-polynomial 
dynamical systems. 



Introduction 

This article is intended to draw the attention of specialists in the fields of 
special functions, combinatorics and nonlinear dynamics on a new class of 
special functions that solve a broad class of nonlinear dynamical systems. 

Let us define an autonomous dynamical system in continuous time as a sys- 
tem of coupled nonlinear ODEs of the form 

3*1 fi(%li ; *^n) (l) 

where the dot denotes the time derivative, the index % runs from 1 to n, the 
dependent variables Xi are real, the functions fa are enough regular to ensure 
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the existence and unicity of solutions, and the boundary conditions are given 
at the initial time. 

For such systems there does not exist a universal structure of the gen- 
eral solution. This is in strong contrast with the particular case of linear 
autonomous dynamical systems for which the general solution of the Cauchy 
problem is given by the exponential of the constant matrix that characterizes 
each of these systems. 

One reason for this lack of universality is the infinite diversity of possible 
functional forms of the functions fi in equations (pQ). Moreover, the nonlin- 
earity of the equations entails a factorial explosion in the coefficients of the 
Taylor series of the solutions: A coefficient of order k involves a sum of k\ 
terms that depend on the specific form of the functions fi. The recursion 
relations between these coefficients are as difficult to solve as the original 
ODEs. Hence, one should not wonder that in most cases it is impossible to 
find a closed form structure for the general k -th order Taylor coefficient of 
such solutions. 

In this article, we show that for the class of quasi-polynomial systems 
defined in the next chapter, a universal form of the general solution to the 
Cauchy problem can be found. Moreover, as has been shown by E.Kerner 
[7], most systems of non-polynomial nonlinear ODEs relevant for physics, 
chemical kinetics and generally for mathematical modeling in natural and 
social sciences can be reshaped in the form of systems of polynomial ODEs. 
The latter, in turn, are a sub-class of the quasi-polynomial systems of ODEs. 
This general property, thus, extends our general solution to most nonlinear 
systems of ODEs defining dynamical systems. 

The quasi-polynomial dynamical systems 

Let us define the sub-class of dynamical systems (pp) for which the functions 
fi are quasi-polynomials, i.e. finite sums of monomials involving powers of 
the dependent variables that can be real numbers. Such monomials we shall 
call quasi-monomials. 

The systems of this class can be cast into a useful standard notation, the 
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so-called quasi-polynomial (QP) representation[U [2]: 

N n 

ii = XiYl Aij I| x k k forz = l,...,n (2) 

3=1 k=l 

where N refers to the number of quasi-monomials in the set of variables Xj 
in the right-hand side, and A and B are constant rectangular matrices re- 
spectively nx N and N xn with entries that are real numbers. The presence 
of the factor X{ in front of the right-hand side of the above equation is es- 
sential and constrains the definition of both matrices. However, it does not 
restrict the generality of the class of systems ([2]) as any polynomial or quasi- 
polynomial can be represented in that form. The presence of the factor Xj 
underlines the importance of the logarithmic time derivative of the functions 
Xi(t) which plays a fundamental role in the present approach. 

Other authors have found independently the same form j2j) and derived 
from it theoretical results in the fields of chemical reactions and ecological 
systems [3 [H [5]. These works were mainly concerned with models describing 
complex networks of interacting entities and they were essentially devoted to 
the study of the stability properties of some of their solutions. More gener- 
ally, the QP class covers almost all the systems of interest for mathematical 
modeling in natural and social sciences. 

A fundamental feature of the above QP representation of a dynamical 
system is its covariance under the group of quasi-monomial transformations. 
These are transformations of the dependent variables defined as follows: 

n 

%i = IJ x k tk for % = 1, . . . , n (3) 

k=l 

where the matrix C is any invertible, constant, real, square n x n matrix. 
In order to avoid unrelevant discussions about the existence of the inverses 
of these transformations and their differentiability we limit our scope to sys- 
tems and to initial conditions such that their solutions remain in the positive 
cone. However, it can be proven that most of the following results continue 
to be valid for systems not fulfilling this restriction. 

It can be easily shown that under transformations ([3|) the system (j2J) 
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becomes: 

N n 

Xi = Aij IJ ^h 3 " forz = l,...,n (4) 

j=l k=l 

with the following rule of transformations for the matrices A and B: 

A = C~ X A (5) 

and 

B = BC (6) 

where the products are matrix products. The identity of form of equations 
((2j) and (jlj) clearly exhibits the covariance of the QP equations under the 
quasi-monomial group of transformations. Let us stress that: 

BA = BA (7) 

in other words, the N x N matrix BA is an invariant of the quasi-monomial 
transformations (jBJ). This means that the whole set of QP-systems is divided 
into equivalence classes labeled by such N x N matrices. The fundamental 
matrix BA is related to the existence of a canonical form in each equivalence 
class as we now show. 

Indeed, under a particular quasi-monomial transformation, any QP sys- 
tem can be brought to a canonical form, the so-called Lotka-Volterra system 
of ODEs [HE11]: 

N 

%i = x i^2 MijXj for i = 1, . . . , N (8) 

3=1 

where A" is the number of monomials of the original QP system ((2j) and the 
square N x N matrix M is equal to the matricial product BA of the two ma- 
trices B and A. There also exists a second canonical form dual of the above 
one, which until now has been much less studied [1] but is also characterized 
by the same matrix BA. 

In a given equivalence class, the passage to the Lotka-Volterra form is 
most easily shown in the particular case of QP systems for which the number 
N of quasi-monomials is equal to the dimension n of the system ([2j), the 
so-called square QP-systems. For the square QP-systems the matrices A and 
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B are of course square matrices. In this case, and if the matrix B is not 
singular, the quasi-monomial transformation ((3]) in which the matrix C is 
chosen equal to B~ x leads to the equation (TjJ with 

A = BA (9) 

and 

B = I (10) 
where / is the identity matrix. Thus, the transformed equation is now: 

N 

x l = x i Y / (BA) ij x j for i = 1, . . . , N (11) 

where we omitted for convenience the tilde accent on the new variables 
Xi in the transformed equation. This is exactly the Lotka-Volterra form an- 
nounced in (JSj) with M = BA, i.e. M is the invariant matrix ((71). 

For the non-square QP systems, i.e. the systems whose matrices A and B 
respectively are nx N and N x n, the Lotka-Volterra form (JSj) is also shown 
to be the canonical form with a matrix M = BA that is N x N [2, 6j. Thus, 
the reduction to the Lotka-Volterra canonical form is a general property of 
the QP dynamical systems. 

The fact that the transformation leading to the canonical form is a dif- 
feomorphism reduces the analysis of the solutions of general QP-systems to 
those of the corresponding Lotka-Volterra systems. This result paves the 
way to the transfer of the whole corpus of knowledge related to the latter 
equations which is pretty extensive. 

We now show that this reduction to the Lotka-Volterra format leads to 
an explicit formula for the general solution of the QP-systems. This solution 
is expressed in terms of a Taylor series whose general coefficient is explicitely 
known. 

General solution of QP-systems 

The Taylor series in time t for the solution Xi{t) of the Lotka-Volterra differ- 
ential system (JSj) is defined as 
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oo ^k 

Zi(*) = £ci(fc)7T for i = 1,...,N (12) 
k=o K[ 

The coefficient q(/c) is calculated by performing the k order time deriva- 
tive of Xi(t) at time t — 0. This is readily done by iterating the t derivative 
while recursively using the Lotka-Volterra system (jEJ). The result is amaz- 
ingly simple: 

N N 

d(k) = xi(0) E • • • E M -i ( M -2 + M ilia ) . . . (M iifc + M hik + ■■■ + M ik _ lik ) 

il=0 ife=0 

x n (0)...x Jfc (0)(13) 

where the Xj(0) are the components of the initial condition. Here again 
we omitted the tilde accent on the variables Xi . 

The considered Lotka-Volterra system results from the transformation 
( l8f9lfT0l) of a QP-system (j2J). Going back to the original variables of this QP- 
system, we get for the order k coefficient Ci{k) of the Taylor series solving 
this system: 

N N 

d{k) = x t (0) E • • • E A ih {Aa 2 + M lll2 ) . . . (Ai ik + M lllk + ... + M tk _ llk ) 

U=0 ifc=0 

n n 

.r ; J()) /; - ... II Mi l! 

ii=i ifc=i 

where % runs from 1 to n and n is the dimension of the original QP-system. 
Thus, the general solution of the QP-system is 

oo 

Xi{t) = J2 C i( k )T, fori = l,...,n (15) 

An observation of the coefficients q(fc) clearly reveals a combinatorial 
structure related to the product 

M ih (M ll2 + M lll2 ) . . . (M iik + M lllk + ... + M tk _ llk ) (16) 

For a one dimensional Lotka-Volterra system this product would reduce 
to k\M k . In more dimensions the number M is replaced by the components 
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of the matrix My. By doubling the number of sums and introducing the 
corresponding Kronecker deltas the previous product becomes 

N N 

X ■ ■ ■ X % + <5uj 2 ) • • • (<% + Kh + ■■■ + K-ih) M hh M j 2 i2 ■ ■ ■ M hi k 

h=0 j k =0 

(17) 

where the tensor 

kh (S ij2 + S hh ) . . . [6 ijh + 5 iljk + ... + S ik _ ljk ) (18) 

is a N-dimensional generalization of the factorial function k\. Its role in the 
structure of the Taylor coefficient is double. It governs the type of contracted 
products between the k tensors M that appear in the coefficient Ci(k), and 
also counts these types. The combinatorial study of this object would cer- 
tainly lead to interesting perspectives. 

The above Taylor series (fl2j) . when it converges, can be considered as 
defining a new class of special functions which presents some analogy with 
the large class of the hypergeometric functions. In analogy with the latter 
class of special functions, the study of other properties of these new functions 
such as asymptotic series expansions and integral representations would be 
of the greatest interest. Indeed, the QP-systems are known to have for cer- 
tain values of their parameters chaotic solutions. The latter have extremely 
complex behaviours in t such as the sensitivity with respect to the variations 
of the initial conditions or strange attractors having fractal geometry in the 
phase-space. These properties should be related in some way to the proper- 
ties of the functions that we generated as solutions of these systems. Such a 
study would certainly contribute to develop a bridge between the theory of 
nonlinear dynamical systems and the theory of special functions. 
Furthermore, we believe that the study of the combinatorial properties of the 
tensor (TT8l) could be of great interest for both combinatorics and dynamical 
system theory. 
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